library(methods)
library(XML)
library(corrplot)
library(scales)
library(lubridate)
library(moveHMM)
library(tidyverse)
set.seed(1337)
load("/home/janique/mesa/R/chosen_people.RData")
padded_ids <- str_pad(chosen_ids, width = 4, pad = "0", side = "left")
summary_df <- data.frame(matrix(ncol = 18, nrow = 0))
col_names <- c("ID", "state_mean1", "state_mean2", "state_mean3",
"sleep_integral", "stay_sleep", "stay_awake", "centre_point",
"amplitude", "accuracy_acti", "accuracy_psg",
"nb_rest_states", "median_limit", "mean_active", "mean_rest",
"var_active", "var_rest", "ri")
colnames(summary_df) <- col_names
for (pad_id in padded_ids) {
print(paste("Processing participants", pad_id))
file_name <- paste0("mesa-sleep-", pad_id, ".csv")
summary_df[nrow(summary_df) + 1,] = fit_model(file_name)
}
## [1] "Processing participants 6003"
## [1] "Processing participants 3543"
## [1] "Processing participants 6155"
## [1] "Processing participants 3904"
## [1] "Processing participants 2133"
## [1] "Processing participants 2049"
## [1] "Did not find a HMM for 3 states"
## [1] "NO SUMMARY STATISTICS ADDED, MODEL DID NOT CONVERGE"
## [1] "Processing participants 1717"
## [1] "Processing participants 4849"
## [1] "Processing participants 6754"
## [1] "Did not find a HMM for 3 states"
## [1] "NO SUMMARY STATISTICS ADDED, MODEL DID NOT CONVERGE"
## [1] "Processing participants 3024"
## [1] "Processing participants 5792"
## [1] "Processing participants 6697"
## [1] "Processing participants 5875"
## [1] "Processing participants 3495"
## [1] "Processing participants 3737"
## [1] "Processing participants 3572"
## [1] "Processing participants 6567"
## [1] "Processing participants 2987"
## [1] "Processing participants 1058"
## [1] "Processing participants 1026"
## [1] "Processing participants 0838"
## [1] "Processing participants 5734"
## [1] "Processing participants 5697"
## [1] "Processing participants 6804"
## [1] "Processing participants 0474"
## [1] "Processing participants 2233"
## [1] "Processing participants 5369"
## [1] "Processing participants 2709"
## [1] "Processing participants 0672"
## [1] "Processing participants 4333"
## [1] "Processing participants 4477"
## [1] "Processing participants 3629"
## [1] "Processing participants 4794"
## [1] "Did not find a HMM for 3 states"
## [1] "NO SUMMARY STATISTICS ADDED, MODEL DID NOT CONVERGE"
## [1] "Processing participants 3674"
## [1] "Did not find a HMM for 3 states"
## [1] "NO SUMMARY STATISTICS ADDED, MODEL DID NOT CONVERGE"
## [1] "Processing participants 1742"
## [1] "Processing participants 4478"
## [1] "Processing participants 2665"
## [1] "Processing participants 4408"
## [1] "Processing participants 4168"
## [1] "Processing participants 6036"
## [1] "Processing participants 4763"
## [1] "Processing participants 6317"
## [1] "Processing participants 5993"
## [1] "Processing participants 5932"
## [1] "Processing participants 6022"
## [1] "Processing participants 4109"
## [1] "Processing participants 6557"
## [1] "Processing participants 6372"
## [1] "Processing participants 2172"
## [1] "Processing participants 5396"
for (i in seq_along(summary_df$nb_rest_states)) {
if (is.na(summary_df$nb_rest_states[i])) {
summary_df[i, ] <- fit_model(paste0("mesa-sleep-", summary_df$ID[i],
".csv"))
}
}
save(summary_df, file = "/home/janique/mesa/R/summary_df20.RData")
# load("/home/janique/mesa/R/summary_df30.RData")
4 did not converge (out of 50) at 20. 2 did not converge at 25. 1 does not converge at 30.
2-state HMM paper that utilises only high quality data has better sleepers.
## [1] "Statistics from 2-state HMM paper"
## # A tibble: 1 × 5
## trouble_sleep wakeup age sleepy back_sleep
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.61 3.12 68.4 1.94 2.15
## [1] "Statistics for all MESA participants"
## # A tibble: 1 × 5
## trouble_sleep wakeup age sleepy back_sleep
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2.06 3.54 69.6 2.04 2.44
## [1] "Statistics for my chosen patients"
## # A tibble: 1 × 5
## trouble_sleep wakeup age sleepy back_sleep
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2.44 3.8 73.4 2.02 2.24
## [1] "Mean accuracy for my HMM vs PSG: 0.78"
## [1] "50 participants were analysed."
4. Older people are less active
##
## Pearson's product-moment correlation
##
## data: rel_person$sleepage5c and rel_person$state_mean3
## t = -2.0513, df = 48, p-value = 0.02286
## alternative hypothesis: true correlation is less than 0
## 95 percent confidence interval:
## -1.00000000 -0.05194707
## sample estimates:
## cor
## -0.2839014
6. Centre of sleep
## # A tibble: 50 × 2
## mesaid ri
## <dbl> <dbl>
## 1 2133 0.13
## 2 5792 0.29
## 3 1717 0.3
## 4 6155 0.36
## 5 2665 0.4
## 6 6557 0.42
## 7 1026 0.46
## 8 4408 0.48
## 9 2709 0.49
## 10 6003 0.49
## 11 3543 0.51
## 12 6036 0.52
## 13 3629 0.53
## 14 6567 0.53
## 15 4478 0.54
## 16 2172 0.56
## 17 3495 0.56
## 18 2987 0.57
## 19 5993 0.57
## 20 6804 0.57
## 21 5369 0.59
## 22 3572 0.6
## 23 4763 0.6
## 24 5697 0.61
## 25 4794 0.63
## 26 3674 0.64
## 27 5932 0.67
## 28 672 0.68
## 29 838 0.69
## 30 4168 0.69
## 31 6754 0.69
## 32 4109 0.7
## 33 2233 0.72
## 34 5734 0.72
## 35 4333 0.73
## 36 474 0.74
## 37 6372 0.74
## 38 4849 0.76
## 39 3024 0.78
## 40 3737 0.78
## 41 3904 0.84
## 42 4477 0.84
## 43 6022 0.84
## 44 2049 0.85
## 45 6317 0.86
## 46 6697 0.86
## 47 1058 0.89
## 48 5396 0.91
## 49 5875 0.96
## 50 1742 NA
7. Level of activity and circadian rhythm
##
## Pearson's product-moment correlation
##
## data: rel_person$mean_active and rel_person$ri
## t = 3.3453, df = 47, p-value = 0.0008113
## alternative hypothesis: true correlation is greater than 0
## 95 percent confidence interval:
## 0.2240296 1.0000000
## sample estimates:
## cor
## 0.4385346
8. Regression
# lm1 <- lm(ri ~ sleepage5c*mean_active + as.factor(gender1), data = rel_person)
# summary(lm1)